'''
Created on Sep 28, 2009

@author: mkiyer
'''
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import sys
import operator

def plot_coverage_profile(utrexons1, utrexons2):
    """
    Make a histogram of normally distributed random numbers and plot the
    analytic PDF over it
    """
    
    x1 = []
    for e in utrexons1:
        if e.max >= 2.0:
            x1.append(e.max)

    x2 = []
    for e in utrexons2:
        if e.max >= 2.0:
            x2.append(e.max)


    fig = plt.figure()
    ax = fig.add_subplot(111)
    nbins=100
    
    # the histogram of the data
    n, bins, patches = ax.hist(x1, bins=nbins, facecolor='green', alpha=0.25)
    n, bins, patches = ax.hist(x2, bins=bins, facecolor='red', alpha=0.25)

    ax.grid(True)
    plt.show()
    return
    
    mu, sigma = 100, 15
    x = mu + sigma * np.random.randn(10000)
    
    fig = plt.figure()
    ax = fig.add_subplot(111)
    
    # the histogram of the data
    n, bins, patches = ax.hist(x, 50, normed=1, facecolor='green', alpha=0.75)
    
    # hist uses np.histogram under the hood to create 'n' and 'bins'.
    # np.histogram returns the bin edges, so there will be 50 probability
    # density values in n, 51 bin edges in bins and 50 patches.  To get
    # everything lined up, we'll compute the bin centers
    bincenters = 0.5*(bins[1:]+bins[:-1])
    # add a 'best fit' line for the normal PDF
    y = mlab.normpdf( bincenters, mu, sigma)
    l = ax.plot(bincenters, y, 'r--', linewidth=1)
    
    ax.set_xlabel('Smarts')
    ax.set_ylabel('Probability')
    #ax.set_title(r'$\mathrm{Histogram\ of\ IQ:}\ \mu=100,\ \sigma=15$')
    ax.set_xlim(40, 160)
    ax.set_ylim(0, 0.03)
    ax.grid(True)
    
    plt.show()
    plt.close()


def plot_exon_profile(utrexons):

    fig = plt.figure()
    ax = fig.add_subplot(111)
    
    x = []
    y = []
    for e in utrexons:
        if e.max > 100.0:
            x.append(math.log(e.end - e.start, 10))
            y.append(math.log(e.mean, 10))
    
    ax.scatter(x,y)
    ax.grid(True)
    plt.xlabel('log10 exon length')
    plt.ylabel('log10 exon coverage')
    plt.show()
    plt.close()
    pass

def read_exon_profile(fhd):
    UTRExon = type('UTRExon', (object,), dict())

    for linenum, line in enumerate(fhd):
        line = line.strip()
        if line is None:
            continue
        if line.startswith('#'):
            continue
        fields = line.split('\t')
        ue = UTRExon()        
        ue.chrom = fields[0]
        ue.start = int(fields[1])
        ue.end = int(fields[2])
        ue.name = fields[3]
        ue.min = float(fields[4])
        ue.max = float(fields[5])
        ue.mean = float(fields[6])
        ue.std = float(fields[7])
        ue.zscore = float(fields[8])
        yield ue    

if __name__ == '__main__':
    
    normal = list(read_exon_profile(open(sys.argv[1])))
    cancer = list(read_exon_profile(open(sys.argv[2])))
    plot_exon_profile(cancer)
    plot_coverage_profile(normal, cancer)

    sorted_cancer = sorted(cancer, key=operator.attrgetter('max'))
    for e in sorted_cancer:
        print e.name, e.chrom, e.start, e.end, 'length', e.end - e.start, e.min, e.max, e.mean, e.std, e.zscore

    
    
    
    pass